import numpy as np


def create(n):
    """
    Create a n-dim Trigonometric function

    :param n: dimensionality
    :return: x0, f, g, G
    """

    x0 = np.array([1.0 / n for i in range(n)], dtype=np.float64)
    idx = np.arange(1, n + 1, dtype=np.float64)

    def fs(x: np.ndarray):
        assert x.ndim == 1
        assert x.shape[0] == n
        sin_x = np.sin(x)
        cos_x = np.cos(x)
        part = np.sum(cos_x)
        fx = (idx - idx * cos_x - sin_x) + (n - part)
        return fx, sin_x, cos_x

    def f(x: np.ndarray):
        return np.sum(fs(x)[0] ** 2)

    def g(x: np.ndarray):
        fx, sin_x, cos_x = fs(x)
        fx = fx * 2
        return fx * (idx * sin_x - cos_x) + (fx.sum() * sin_x)

    def h(x: np.ndarray):
        fx, sin_x, cos_x = fs(x)
        isc = idx * sin_x - cos_x
        d_isc = idx * cos_x + sin_x
        # some NumPy-broadcasting magic
        dfx = np.diag(isc) + sin_x
        s_dfx = dfx.sum(axis=0)
        h1 = np.diag(fx.sum() * cos_x + fx * d_isc)
        h2 = dfx * isc[:, np.newaxis]
        h3 = sin_x[:, np.newaxis] @ s_dfx[np.newaxis, :]
        return (h1 + h2 + h3) * 2

    return x0, f, g, h
